﻿Imports System.IO
Imports System.Windows.Forms.DataVisualization.Charting
Imports System.Math
Module DPotModule
    Public DataStr() As String
    Public CmputRes As String
    Dim DataStream As String
    Public PotData() As DPotClass
    Public PotData2() As DPotClass
    Public PotData3() As DPotClass
    Public PotData4() As DPotClass
    Dim PotNum As Integer
    Public a, mf, L0, B0 As Double
    Public mb, e2, e_2, W, n2, t, N, M, M0 As Double
    Public mA1(5) As Double
    Public mA2(6) As Double
    Public Pot1() As pot
    Structure pot
        Dim x, y As Double
        Dim id As String
    End Structure
    Sub OpenFile()
        Dim OPEN As New OpenFileDialog
        OPEN.Filter = "TEXT|*.txt"
        OPEN.FilterIndex = 0
        If OPEN.ShowDialog = DialogResult.OK Then
            Try
                Dim strm As New System.IO.StreamReader(OPEN.FileName)
                DataStream = strm.ReadToEnd
                DataStream = DataStream.Replace(vbCrLf, ",")
                DataStream = DataStream.Replace(vbLf, ",")
                DataStream = DataStream.Trim(",")
                DataStr = DataStream.Split(",")
                PotNum = (DataStr.Length - 9) / 4
                ReDim PotData(PotNum - 1)
                a = DataStr(1)
                mf = DataStr(3)
                L0 = CRad(DataStr(5))
                B0 = CRad(DataStr(7))
                For i = 0 To PotNum - 1
                    PotData(i) = New DPotClass
                    PotData(i).PotName = DataStr(9 + i * 4)
                    PotData(i).B = CRad(DataStr(9 + i * 4 + 1))
                    PotData(i).L = CRad(DataStr(9 + i * 4 + 2))
                    PotData(i).H = DataStr(9 + i * 4 + 3)
                    PotData(i).flag = 1
                Next
            Catch ex As Exception
                MsgBox("未导入数据"， 64， "警告")
                DataStr = Nothing
                Exit Sub
            End Try
        Else
            MsgBox("未选择文件"， 64， "警告")
            Exit Sub
        End If
    End Sub
    Sub SaveFile()
        Dim save As New SaveFileDialog
        save.Filter = "TEXT|*.txt"
        save.FilterIndex = 0
        If save.ShowDialog = DialogResult.OK Then
            Dim strm As New System.IO.StreamWriter(save.OpenFile)
            strm.Write(CmputRes)
            strm.Close()
        Else
            MsgBox("未保存成功"， 64， "警告")
        End If
    End Sub
    Sub bmp(ByVal chart1 As Chart)
        Dim bmp As New SaveFileDialog
        bmp.Filter = "BMP|*.bmp"
        If bmp.ShowDialog = DialogResult.OK Then
            chart1.SaveImage(bmp.FileName, System.Drawing.Imaging.ImageFormat.Bmp)
        Else
            MsgBox("未保存成功"， 64， "警告")
        End If
    End Sub
    Function GetDataGrid()
        Dim Table As New DataTable
        For i = 0 To 3
            Table.Columns.Add()
        Next
        Table.Columns(0).ColumnName = "点名"
        Table.Columns(1).ColumnName = "纬度B"
        Table.Columns(2).ColumnName = "经度L"
        Table.Columns(3).ColumnName = "椭球高（H）"
        For i = 0 To PotNum - 1
            Table.Rows.Add()
            Table(i)(0) = PotData(i).PotName
            Table(i)(1) = RadToDmsStr(PotData(i).B)
            Table(i)(2) = RadToDmsStr(PotData(i).L)
            Table(i)(3) = PotData(i).H
        Next
        Return Table
    End Function
    Sub CalAll()
        GetElliElement()
        For i = 0 To PotNum - 1
            PotData(i).BLHToXYZ()
        Next
        GetXYZ()
        For i = 0 To PotNum - 1
            PotData2(i).XYZToBLH()
        Next

        GetBL()
        For i = 0 To PotNum - 1
            PotData3(i).BLToXY()
        Next
        GetXY()
        For i = 0 To PotNum - 1
            PotData4(i).XYToBL()
        Next
        GetPot1()
    End Sub
    Sub GetPot1()
        ReDim Pot1(PotNum - 1)
        For i = 0 To PotNum - 1
            Pot1(i).x = PotData(i).X
            Pot1(i).y = PotData(i).Y
        Next
    End Sub
    Sub GetElliElement()
        mb = a - a / mf
        e2 = (Pow(a, 2) - Pow(mb, 2)) / Pow(a, 2)
        e_2 = e2 / (1 - e2)
        W = Sqrt(1 - e2 * Pow(Sin(B0), 2))
        n2 = e_2 * Pow(Cos(B0), 2)
        t = Tan(B0)
        N = a / W
        M = a * (1 - e2) / Pow(W, 3)
        M0 = a * (1 - e2)
    End Sub
    Sub GetXYZ()
        ReDim PotData2(PotNum - 1)
        For i = 0 To PotNum - 1
            PotData2(i) = New DPotClass
            PotData2(i).PotName = PotData(i).PotName
            PotData2(i).X = PotData(i).X
            PotData2(i).Y = PotData(i).Y
            PotData2(i).Z = PotData(i).Z
            PotData2(i).X += 2018 : PotData2(i).Y += 2018 : PotData2(i).Z += 2018
            PotData2(i).flag = 2
        Next
    End Sub
    Sub GetXY()
        ReDim PotData4(PotNum - 1)
        For i = 0 To PotNum - 1
            PotData4(i) = New DPotClass
            PotData4(i).PotName = PotData3(i).PotName
            PotData4(i).X = PotData3(i).X
            PotData4(i).Y = PotData3(i).Y
            PotData4(i).X += 1000 : PotData4(i).Y += 1000
            PotData4(i).flag = 4
        Next
    End Sub
    Sub GetBL()
        ReDim PotData3(PotNum - 1)
        For i = 0 To PotNum - 1
            PotData3(i) = New DPotClass
            PotData3(i).PotName = PotData(i).PotName
            PotData3(i).B = PotData(i).B
            PotData3(i).L = PotData(i).L
            PotData3(i).flag = 3
        Next
    End Sub
    Function GetArray1(pm As Double)
        Dim mAa(5) As Double
        Dim Ac, Bc, Cc, Dc, Ec, Fc As Double
        Ac = 1 + 3 / 4 * e2 + 45 / 64 * Pow(e2, 2) + 175 / 256 * Pow(e2, 3) + 11025 / 16384 * Pow(e2, 4) + 43659 / 65536 * Pow(e2, 5)
        Bc = 3 / 4 * e2 + 15 / 16 * Pow(e2, 2) + 525 / 512 * Pow(e2, 3) + 2205 / 2048 * Pow(e2, 4) + 72765 / 65536 * Pow(e2, 5)
        Cc = 15 / 64 * Pow(e2, 2) + 105 / 256 * Pow(e2, 3) + 2205 / 4096 * Pow(e2, 4) + 10395 / 16384 * Pow(e2, 5)
        Dc = 35 / 512 * Pow(e2, 3) + 315 / 2048 * Pow(e2, 4) + 31185 / 131072 * Pow(e2, 5)
        Ec = 315 / 16384 * Pow(e2, 4) + 3465 / 65536 * Pow(e2, 5)
        Fc = 693 / 131072 * Pow(e2, 5)
        mAa(0) = Ac * pm
        mAa(1) = -1 / 2 * Bc * pm
        mAa(2) = 1 / 4 * Cc * pm
        mAa(3) = -1 / 6 * Dc * pm
        mAa(4) = 1 / 8 * Ec * pm
        mAa(5) = -1 / 10 * Fc * pm
        Return mAa
    End Function
    Function GetArray2(pB As Double)
        Dim mAa(6) As Double
        Dim mt As Double = Getf(pB)(3)
        Dim pmn As Double = Getf(pB)(0)
        Dim mn As Double = Getf(pB)(1)
        mA1 = GetArray1(M0)
        mAa(0) = mA1(0) * pB + mA1(1) * Sin(2 * pB) + mA1(2) * Sin(4 * pB) + mA1(3) * Sin(6 * pB) + mA1(4) * Sin(8 * pB) + mA1(5) * Sin(10 * pB)
        mAa(1) = pmn * Cos(pB)
        mAa(2) = 1 / 2 * pmn * Pow(Cos(pB), 2) * mt
        mAa(3) = 1 / 6 * pmn * Pow(Cos(pB), 3) * (1 - mt ^ 2 + mn ^ 2)
        mAa(4) = 1 / 24 * pmn * Pow(Cos(pB), 4) * (5 - mt ^ 2 + 9 * mn ^ 2 + 4 * mn ^ 4) * mt
        mAa(5) = 1 / 120 * pmn * Pow(Cos(pB), 5) * (5 - 18 * mt ^ 2 + mt ^ 4 + 14 * mn ^ 2 - 58 * mn ^ 2 * mt ^ 2)
        mAa(6) = 1 / 720 * pmn * Pow(Cos(pB), 6) * (61 - 58 * mt ^ 2 + mt ^ 4 + 270 * mn ^ 2 - 330 * mn ^ 2 * mt ^ 2) * mt
        Return mAa
    End Function
    Function GetBf(pX As Double)
        Dim Bf As Double
        Dim mA1 = GetArray1(M0)
        Dim pB0 As Double
        Bf = pX / mA1(0)
        Dim Deta As Double
        Do
            pB0 = Bf
            Deta = mA1(1) * Sin(2 * pB0) + mA1(2) * Sin(4 * pB0) + mA1(3) * Sin(6 * pB0) + mA1(4) * Sin(8 * pB0) + mA1(5) * Sin(10 * pB0)
            Bf = (pX - Deta) / mA1(0)
        Loop While Abs(Bf - pB0) > 0.00000001
        Return Bf
    End Function
    Function Getf(pBf As Double)
        Dim ma(4) As Double
        Dim Wf, Nf, nf_, Mf, tf As Double
        Wf = Sqrt(1 - e2 * Pow(Sin(pBf), 2))
        nf_ = Sqrt(e_2 * Pow(Cos(pBf), 2))
        tf = Tan(pBf)
        Nf = a / Wf
        Mf = a * (1 - e2) / Pow(Wf, 3)
        ma(0) = Nf
        ma(1) = nf_
        ma(2) = Mf
        ma(3) = tf
        ma(4) = pBf
        Return ma
    End Function
    Function GetArray3(ma() As Double)
        Dim ArrayB(6) As Double
        Dim Nf, nf_, Mf, tf, Bf As Double
        Nf = ma(0)
        nf_ = ma(1)
        Mf = ma(2)
        tf = ma(3)
        Bf = ma(4)
        ArrayB(0) = Bf
        ArrayB(1) = 1 / Nf / Cos(Bf)
        ArrayB(2) = -tf / 2 / Nf / Mf
        ArrayB(3) = -(1 + 2 * tf ^ 2 + nf_ ^ 2) * ArrayB(1) / (6 * Nf ^ 2)
        ArrayB(4) = -(5 + 3 * tf ^ 2 + nf_ ^ 2 - 9 * nf_ ^ 2 * tf ^ 2) * ArrayB(2) / (12 * Nf ^ 2)
        ArrayB(5) = -(5 + 28 * tf ^ 2 + 24 * tf ^ 4 + 6 * nf_ ^ 2 + 8 * nf_ ^ 2 * tf ^ 2) * ArrayB(1) / (120 * Nf ^ 4)
        ArrayB(6) = (61 + 90 * tf ^ 2 + 45 * tf ^ 4) * ArrayB(2) / (360 * Nf ^ 4)
        Return ArrayB
    End Function
    Function CRad(ByVal Dms As Double)
        Dim ValueSign As Integer
        Dim DD, MM As Integer
        Dim SS As Double
        ValueSign = Sign(Dms)
        Dms = ValueSign * Dms
        DD = Int(Dms)
        MM = Int((Dms - DD) * 100)
        SS = (((Dms - DD) * 100 - MM) * 100)
        Return (DD + MM / 60 + SS / 3600) * ValueSign / 180 * PI
    End Function
    Function RadToDmsStr(ByVal Rad As Double)
        Dim DMS As Double
        Dim DD, MM As Integer
        Dim SS As Double
        Dim ValueSign As Integer
        ValueSign = Sign(Rad)
        Rad = ValueSign * Rad
        DMS = Rad / PI * 180
        DD = Int(DMS)
        MM = Int((DMS - DD) * 60)
        SS = ((DMS - DD) * 60 - MM) * 60
        If SS > 60 Then
            SS -= 60
            MM += 1
        End If
        If MM > 60 Then
            MM -= 60
            DD += 1
        End If
        If ValueSign = -1 Then
            Return String.Format("-{0}°{1:00}′{2:00.0000}″", DD, MM, SS)
        Else
            Return String.Format("{0}°{1:00}′{2:00.0000}″", DD, MM, SS)
        End If
    End Function
    Sub SaveBMP(ByVal Chart1 As Chart) 'BMP图像保存过程
        Dim SFD As New SaveFileDialog
        SFD.Filter = "BMP格式|*.bmp"
        SFD.FilterIndex = 0
        If SFD.ShowDialog = DialogResult.OK Then
            Chart1.SaveImage(SFD.FileName, System.Drawing.Imaging.ImageFormat.Bmp)
            MsgBox("保存完成", 64, "提示")
        End If
    End Sub
    Sub DXF_start()
        PrintLine(1, 0)
        PrintLine(1, "SECTION")
        PrintLine(1, 2)
        PrintLine(1, "ENTITIES")
    End Sub
    Sub DXF_end()
        PrintLine(1, 0)
        PrintLine(1, "ENDSEC")
        PrintLine(1, 0)
        PrintLine(1, "EOF")
    End Sub
    Sub DXF_Point(ByVal x As Double, ByVal y As Double)
        PrintLine(1, 0)
        PrintLine(1, "POINT")
        PrintLine(1, 8)
        PrintLine(1, "points")
        PrintLine(1, 62)
        PrintLine(1, "1")
        PrintLine(1, 10)
        PrintLine(1, x)
        PrintLine(1, 20)
        PrintLine(1, y)
        PrintLine(1, 30)
        PrintLine(1, "0.0")
    End Sub
    Sub DXF_Line(ByVal x1 As Double, ByVal y1 As Double, ByVal x2 As Double, ByVal y2 As Double)
        PrintLine(1, 0)
        PrintLine(1, "LINE")
        PrintLine(1, 8)
        PrintLine(1, "LINEs")
        PrintLine(1, 62)
        PrintLine(1, "1")
        PrintLine(1, 10)
        PrintLine(1, x1)
        PrintLine(1, 20)
        PrintLine(1, y1)
        PrintLine(1, 11)
        PrintLine(1, x2)
        PrintLine(1, 21)
        PrintLine(1, y2)
    End Sub
    Sub DXF_Text(ByVal x As Double, ByVal y As Double, ByVal Name As String)
        PrintLine(1, 0)
        PrintLine(1, "TEXT")
        PrintLine(1, 8)
        PrintLine(1, "POTNAME")
        PrintLine(1, 62)
        PrintLine(1, 4)
        PrintLine(1, 100)
        PrintLine(1, "ACDBTEXT")
        PrintLine(1, 10)
        PrintLine(1, x)
        PrintLine(1, 20)
        PrintLine(1, y)
        PrintLine(1, 30)
        PrintLine(1, "0.0")
        PrintLine(1, 40)
        PrintLine(1, 2)
        PrintLine(1, 1)
        PrintLine(1, Name)
        PrintLine(1, 7)
        PrintLine(1, "宋体")
    End Sub
    Sub DrawDXF()
        Dim SFD As New SaveFileDialog
        SFD.Filter = "DXF文件(*.dxf)|*.dxf"
        SFD.FilterIndex = 0
        If SFD.ShowDialog = DialogResult.OK Then
            FileOpen(1, SFD.FileName, OpenMode.Output)
            DXF_start()

            For i = 0 To Pot1.Length - 1
                DXF_Point(Pot1(i).x, Pot1(i).y)
                DXF_Text(Pot1(i).x, Pot1(i).y, Pot1(i).id)
            Next

            DXF_end()
            FileClose(1)
            MsgBox("导出完成", 64, "提示")
        End If
    End Sub
    Sub Draw(ByVal Chart1 As Chart)
        Dim Xmax, Xmin, Ymax, Ymin As Double
        Dim X(Pot1.Length - 1), Y(Pot1.Length - 1) As Double
        For i = 0 To Pot1.Length - 1
            X(i) = Pot1(i).x
            Y(i) = Pot1(i).y
        Next
        Xmax = X.Max : Xmin = X.Min
        Ymax = Y.Max : Ymin = Y.Min

        Chart1.ChartAreas.Clear()
        Chart1.Series.Clear()
        Chart1.Titles.Clear()

        Dim Tit As New Title("散点图")
        Tit.Font = New Font("黑体", 12.5, FontStyle.Regular)
        Chart1.Titles.Add(Tit)

        Dim Area1 As New ChartArea
        Area1.AxisX.MajorTickMark.Enabled = False
        Area1.AxisY.MajorTickMark.Enabled = False
        Area1.AxisX.Title = ""
        Area1.AxisY.Title = ""
        Area1.AxisX.Minimum = (Xmin \ 10 - 1) * 10
        Area1.AxisY.Minimum = (Ymin \ 10 - 1) * 10
        Area1.AxisX.Maximum = (Xmax \ 10 + 1) * 10
        Area1.AxisY.Maximum = (Ymax \ 10 + 1) * 10
        Chart1.ChartAreas.Add(Area1)


        Dim Series1(Pot1.Length - 1) As Series
        For i = 0 To Pot1.Length - 1
            Series1(i) = New Series
            Series1(i).ChartType = SeriesChartType.Point
            Series1(i).MarkerSize = 10
            Series1(i).MarkerStyle = MarkerStyle.Circle
            Series1(i).MarkerColor = Color.Black
            Series1(i).IsVisibleInLegend = False
            Series1(i).Points.AddXY(Pot1(i).x, Pot1(i).y)
            Series1(i).Label = Pot1(i).id
            Chart1.Series.Add(Series1(i))
        Next
    End Sub
    Function Report()
        Dim RText As String
        RText = "1.椭球相关参数计算" & vbCrLf
        RText += String.Format("f:{0:f6}{1}", 1 / mf, vbCrLf)
        RText += String.Format("e^2:{0:f6}{1}", e2, vbCrLf)
        RText += String.Format("e^'2:{0:f6}{1}", e_2, vbCrLf)
        RText += String.Format("W:{0:f6}{1}", W, vbCrLf)
        RText += String.Format("n^2:{0:f6}{1}", n2, vbCrLf)
        RText += String.Format("t:{0:f6}{1}", t, vbCrLf)
        RText += String.Format("N:{0:f6}{1}", N, vbCrLf)
        RText += String.Format("M:{0:f6}{1}", M, vbCrLf)
        RText += String.Format("M0:{0:f6}{1}", M0, vbCrLf)
        RText += vbCrLf
        RText += "2.大地坐标（BLH）转换为空间坐标（XYZ）" & vbCrLf
        RText += "-------------------------" & vbCrLf
        RText += String.Format("{0,-5}{1}{2,-20}{3}{4,-20}{5}{6,-10}{7}{8,-20:f4}{9}{10,-20:f4}{11}{12,-20:f4}{13}", "点名", vbTab, "B", vbTab, "L", vbTab, "H", vbTab, "X", vbTab, "Y", vbTab, "Z", vbCrLf)
        For i = 0 To PotNum - 1
            With PotData(i)
                RText += String.Format("{0,-5}{1}{2,-20}{3}{4,-20}{5}{6,-10:f4}{7}{8,-20:f4}{9}{10,-20:f4}{11}{12,-20:f4}{13}", .PotName, vbTab, RadToDmsStr(.B), vbTab, RadToDmsStr(.L), vbTab, .H, vbTab, .X, vbTab, .Y, vbTab, .Z, vbCrLf)
            End With
        Next
        RText += vbCrLf
        RText += "3.空间坐标（XYZ）转换为大地坐标（BLH）" & vbCrLf
        RText += "-------------------------" & vbCrLf
        RText += String.Format("{0,-5}{1}{2,-20:f4}{3}{4,-20:f4}{5}{6,-10}{7}{8,-20}{9}{10,-20}{11}{12,-20:f4}{13}", "点名", vbTab, "X", vbTab, "Y", vbTab, "Z", vbTab, "B", vbTab, "L", vbTab, "H", vbCrLf)
        For i = 0 To PotNum - 1
            With PotData2(i)
                RText += String.Format("{0,-5}{1}{2,-20:f4}{3}{4,-20:f4}{5}{6,-10:f4}{7}{8,-20}{9}{10,-20}{11}{12,-20:f4}{13}", .PotName, vbTab, .X, vbTab, .Y, vbTab, .Z, vbTab, RadToDmsStr(.B), vbTab, RadToDmsStr(.L), vbTab, .H, vbCrLf)
            End With
        Next
        RText += vbCrLf
        RText += "4.1.高斯正算（BL-->xy）" & vbCrLf
        RText += "-------------------------" & vbCrLf
        RText += String.Format("{0,-5}{1}{2,-20:f4}{3}{4,-20:f4}{5}{6,-20:f4}{7}{8,-20:f4}{9}", "点名", vbTab, "B", vbTab, "L", vbTab, "X", vbTab, "Y", vbCrLf)
        For i = 0 To PotNum - 1
            With PotData3(i)
                RText += String.Format("{0,-5}{1}{2,-20:f4}{3}{4,-20:f4}{5}{6,-20:f4}{7}{8,-20:f4}{9}", .PotName, vbTab, RadToDmsStr(.B), vbTab, RadToDmsStr(.L), vbTab, .X, vbTab, .Y, vbCrLf)
            End With
        Next
        RText += vbCrLf
        RText += "公式（7）计算结果：" & vbCrLf
        RText += "-------------------------" & vbCrLf
        RText += String.Format("α：{0}{1}", mA1(0), vbCrLf)
        RText += String.Format("β：{0}{1}", mA1(1), vbCrLf)
        RText += String.Format("γ：{0}{1}", mA1(2), vbCrLf)
        RText += String.Format("δ：{0}{1}", mA1(3), vbCrLf)
        RText += String.Format("ε：{0}{1}", mA1(4), vbCrLf)
        RText += String.Format("ζ：{0}{1}", mA1(5), vbCrLf)
        RText += vbCrLf
        RText += "公式（10）计算结果：" & vbCrLf
        RText += "-------------------------" & vbCrLf
        For i = 0 To PotNum - 1
            Dim mpB As Double
            mpB = PotData3(i).B
            mA2 = GetArray2(mpB)
            RText += String.Format("{0,-5}{1}a0:{2,-20:f6}{3}a1:{4,-20:f6}{5}a2:{6,-20:f6}{7}a3:{8,-20:f6}{9}a4:{10,-20:f6}{11}a5:{12,-20:f6}{13}a6:{14,-20:f6}{15}", PotData3(i).PotName, vbTab, mA2(0), vbTab, mA2(1), vbTab, mA2(2), vbTab, mA2(3), vbTab, mA2(4), vbTab, mA2(5), vbTab, mA2(6), vbCrLf)
        Next

        RText += vbCrLf
        RText += "5.1.高斯反算（xy-->BL）" & vbCrLf
        RText += "-------------------------" & vbCrLf
        RText += String.Format("{0}{1}{2}{3}{4}{5}{6}{7}{8}{9}", "点名", vbTab, "X", vbTab, "Y", vbTab, "B", vbTab, "L", vbCrLf)
        For i = 0 To PotNum - 1
            With PotData4(i)
                RText += String.Format("{0,-5}{1}{2,-20:f4}{3}{4,-20:f4}{5}{6,-20:f4}{7}{8,-20:f4}{9}", .PotName, vbTab, .X, vbTab, .Y, vbTab, RadToDmsStr(.B), vbTab, RadToDmsStr(.L), vbCrLf)
            End With
        Next
        RText += vbCrLf
        RText += "公式（13）计算结果：" & vbCrLf
        RText += "-------------------------" & vbCrLf
        For i = 0 To PotNum - 1
            Dim mpB As Double = GetBf(PotData4(i).X)
            Dim mab() As Double
            mab = GetArray3(Getf(mpB))
            RText += String.Format("{0,-5}{1}a0:{2,-20:f6}{3}a1:{4,-20}{5}a2:{6,-20}{7}a3:{8,-20}{9}a4:{10,-20}{11}a5:{12,-20}{13}a6:{14,-20}{15}", PotData3(i).PotName, vbTab, mab(0), vbTab, mab(1), vbTab, mab(2), vbTab, mab(3), vbTab, mab(4), vbTab, mab(5), vbTab, mab(6), vbCrLf)
        Next

        RText += vbCrLf
        Return RText
    End Function
End Module
